Настройка вывода данных:

knitr::opts_chunk$set(echo = TRUE, warning = FALSE)
options(width = 100) # ширина текстового вывода
options(digits = 6) # число знаков после запятой в выводе

Подключение необходимых библиотек:

library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(ggrepel)
## Loading required package: ggplot2
library(readr)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.0     v stringr 1.4.0
## v purrr   0.3.2     v forcats 0.4.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date()        masks base::date()
## x dplyr::filter()          masks stats::filter()
## x lubridate::intersect()   masks base::intersect()
## x dplyr::lag()             masks stats::lag()
## x lubridate::setdiff()     masks base::setdiff()
## x lubridate::union()       masks base::union()
library(ggplot2)
library(ggfortify)
library(modelr)
library(broom)
## 
## Attaching package: 'broom'
## The following object is masked from 'package:modelr':
## 
##     bootstrap
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa
library(forcats)
library(tidyr)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
library(dplyr)
library(memisc)
## Loading required package: lattice
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'memisc'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:dplyr':
## 
##     collect, recode, rename, syms
## The following object is masked from 'package:purrr':
## 
##     %@%
## The following object is masked from 'package:ggplot2':
## 
##     syms
## The following object is masked from 'package:lubridate':
## 
##     is.interval
## The following objects are masked from 'package:stats':
## 
##     contr.sum, contr.treatment, contrasts
## The following object is masked from 'package:base':
## 
##     as.array
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Описание переменных:

1) Подготовка исходных данных:

Загрузка исходных данных

bike <- read_csv ("C:/Users/Юра/Desktop/KP-NIS-R/ПРОЕкт/Проек/bike.csv")
## Parsed with column specification:
## cols(
##   datetime = col_datetime(format = ""),
##   season = col_double(),
##   holiday = col_double(),
##   workingday = col_double(),
##   weather = col_double(),
##   temp = col_double(),
##   atemp = col_double(),
##   humidity = col_double(),
##   windspeed = col_double(),
##   casual = col_double(),
##   registered = col_double(),
##   count = col_double()
## )

Вывод: была загружена база данных по арендам велосипедов. Она содержит в себе 10886 наблюдений, которые представляют собой акты аренды велосипедов с предшествующими ими внешними условями, такими как погодные условия, сезон, календарный день и т.д - над ней будет проводиться анализ и поиск подходящей модели прогнозирования. Нам необходимо подготовить данные к анализу и разбить базу данных на обучающую и тестовую выборки случайным образом, и именно к тестовой выборке будет применяться выявленная на обучающей выборке лучшая модель для прогнозирования значений актов аренды велосипедов.

Подготовка данных к анализу

Удалим переменные “registered” и “casual” в нашей выборке, потому что они являются лишь слагаемыми нашей целей переменной “count”, поэтому не несут в себе большой значимости и могут изменить модели корреляции по причиние их сильной зависимости от целевой переменной:

bike <- bike[,-(10:11)]

Удалим времена года и выделим из дат месяца в обучающей и тестовой выборках, а также время, чтобы повысить точность нашего анализа, который может быть основан на месяцах, а не просто временах года:

## Удаление времён года
bike <- bike[,-2]

## Создание отдельных столбцов с временем и датой в выборках
bike$month <- as.Date(bike$datetime) 
bike$time <- format(as.POSIXct(bike$datetime) ,format = "%H:%M:%S")

## Удаление объединённого столбца в выборках
bike <- bike[,-1]

## Выделение месяцев из дат
bike$month <- format(bike$month, "%m")

## Выделение часов из времени
bike$time <- substr(bike$time, start = 1, stop = 2)

Изменение типа данных даты и времени

bike$month <- as.numeric(as.character(bike$month))
bike$time <- as.numeric(as.character(bike$time))

Разраничение времён суток

bike <- mutate (bike,
                daytime = ifelse(time >= 0, ifelse(time<=6, 1,ifelse(time<=12, 2,ifelse (time <=17,3, 4))), 0))

## Удаление ненужного столбца с часами
bike <- bike[,-10]

Описание новых переменных

  • holiday: выходной (1 - да, 0 - нет);
  • workingday: рабочий день (1 - да, 0 -нет);
  • weather: качественная характеристика погодных условий:
    1: Ясно, малооблачно, немного облачно; 2: Туманно и облачно, малооблачно, переменнооблачно; 3: Лёгкий снегопад, лёгкий дождь с громом и перистыми облакми, лёгкий дождь с перистыми облаками; 4: Ливень, град, гром, туман, снег;
  • temp: температура в цельсиях;
  • atemp: температура в цельсиях (по ощущениям);
  • humidity: относительная влажность;
  • windspeed: скорость ветра;
  • month: месяц года (от 1 до 12);
  • daytime: время суток: 1: ночь; 2: утро; 3: день; 4: вечер;
  • count - целевая переменная: общее число актов аренды;

Получение обучающей и тестовой выборки

sample <- floor(0.8*nrow(bike))

set.seed(10886)

sample_ind <- sample(seq_len(nrow(bike)),size = sample)

bike_train <- bike[sample_ind,]
bike_test <- bike[-sample_ind,]

Прежде чем начать моделирование, необходимо выяснить: как сильно на корреляцию влияют удаленные столбцы registered и casual.Так как целевая переменная count = registered+casual, очевидно, что значение двух нецелевых переменных влияют напрямую на целевую. Задание по прогнозированию: прогноз переменной count. Поэтому, при анализе и моделировании необходимо исключить registered и casual. Однако, чтобы не потерять скрытые зависимости, которые несут в себе удаленные столбцы, проведем анализ этих столбцов на зависимости. В лучшем исходе registered и casual должны показывать одинаковые зависимости и тренды, что позволи потсроить модель без них и без ухудшения самой модели.

2) Разведоточный анализ:

Изучение взаимосвязей факторов с числом зарегистрированных/незарегистрированных клиентов для выявления новых закономерностей для подбора переменных в моделях

Загрузим изначальную выборку с переменными “registered” и “casual”. Сделаем отношение в 0.75, которое означает, что если в одном наблюдении более 75% арендаторов велосипедов были зарегистрированными, то считается, что во время этого наблюдения были только зарегистрированные кленты, если менее - то только незарегистрированные. Такая большая доля выбрана потому, что число незарегистрированных пользователей почти всегда незначительна.

bike_orig <- read_csv ("C:/Users/Юра/Desktop/KP-NIS-R/ПРОЕкт/Проек/bike.csv")
## Parsed with column specification:
## cols(
##   datetime = col_datetime(format = ""),
##   season = col_double(),
##   holiday = col_double(),
##   workingday = col_double(),
##   weather = col_double(),
##   temp = col_double(),
##   atemp = col_double(),
##   humidity = col_double(),
##   windspeed = col_double(),
##   casual = col_double(),
##   registered = col_double(),
##   count = col_double()
## )
bike_orig <- mutate(bike_orig, relation = ifelse(registered/count >= 0.75, "registered", "casual"))
graph1 <- bike_orig %>%

  ggplot(aes(count, temp, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title = "1", colour = "Тип клиента", x = "Число актов аренды", y = "t воздуха")

graph2 <- bike_orig %>%

  ggplot(aes(count, atemp, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title = "2", colour = "Тип клиента", x = "Число актов аренды", y = "t воздуха по ощущениям")

graph3 <- bike_orig %>%

  ggplot(aes(count, weather, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title = "3", colour = "Тип клиента", x = "Число актов аренды", y = "Cостояние погоды")

graph4 <- bike_orig %>%

  ggplot(aes(count, season, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title = "4", colour = "Тип клиента", x = "Число актов аренды", y = "Сезон")

grid.arrange(graph1, graph2, graph3, graph4, ncol=2)

Вывод:

1 график: можно заметить, что при высокой температуре воздуха существует тенденция для незарегистрированных пользователей для обращения к аренде велосипедов. Для зарегистрированных пользователей температура не так важна - их число актов аренды распределено равномерно;

2 график: закономерности практически идентичные графику №1;

3 график: можно заявить, что во время самой плохой степени погоды берётся крайне мало велосипедов, но вся их аренда производится только зарегистрированными пользователями. В остальных случаях, есть как зарегистрированные, так незарегистрированные пользователи (зарегистрированных по-прежнему больше);

4 график: можно заметить, что во время зимнего сезона арендуют велосипеды практически только зарегистрированные пользователи. Больше всего незарегистрированных пользователей арендуют велосипеды весной и летом;

graph5 <- bike_orig %>%

  ggplot(aes(count, workingday, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs (title = "5", colour = "Тип клиента", x = "Число актов аренды", y = "Будние дни")

graph6 <- bike_orig %>%

  ggplot(aes(count, holiday, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title = "6", colour = "Тип клиента", x = "Число актов аренды", y = "Праздничные дни")

graph7 <- bike_orig %>%

  ggplot(aes(count, windspeed, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title = "7", colour = "Тип клиента", x = "Число актов аренды", y = "Скорость ветра")

graph8 <- bike_orig %>%

  ggplot(aes(count, humidity, color = relation)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title = "8", colour = "Тип клиента", x = "Число актов аренды", y = "Влажность")

grid.arrange(graph5, graph6, graph7, graph8, ncol=2)

Вывод:

5 график: можно заметить, что в выходные дни есть тенденция для незарегистрированных пользователей прибегать к аренде велосипедов, в то время как в буние дни берут велосипеды практически только зарегистрированные пользователи;

6 график: можно заметить, что в праздничные есть тенденция для незарегистрированных пользователей прибегать к аренде велосипедов, в то время как в обычные дни берут велосипеды практически только зарегистрированные пользователи;

7 график: нет практически никакой зависиомсти между зарегистрированностью пользователя и скоростью ветра. Все типы клиентов распределны случайным образом;

8 график: можно заметить, что при низкой влажности больше незарегистрированных пользователей берёт велосипеды в аренду. Тем не менее, эта зависимость выражена слабо - можно сказать, что влажность тоже не влияет на тенднецию типр=ов клиентов к аренде велосипедов;

Вывод: анализ зависимостей показал значительные несовпадения зависимостей переменных с рабочими днями, сезоном (месяцы) и влажность. Используя эти переменные при моделировании стоит быть аккуратными, т.к. их применение может исказить реальную ситуацию, что даст большую ошибку при прогнозирования.

Статистический анализ выборок

summary (bike_train)
##     holiday         workingday       weather          temp           atemp          humidity    
##  Min.   :0.0000   Min.   :0.000   Min.   :1.00   Min.   : 0.82   Min.   : 0.76   Min.   :  0.0  
##  1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:1.00   1st Qu.:13.94   1st Qu.:16.66   1st Qu.: 47.0  
##  Median :0.0000   Median :1.000   Median :1.00   Median :20.50   Median :24.24   Median : 62.0  
##  Mean   :0.0286   Mean   :0.679   Mean   :1.42   Mean   :20.26   Mean   :23.68   Mean   : 61.9  
##  3rd Qu.:0.0000   3rd Qu.:1.000   3rd Qu.:2.00   3rd Qu.:26.24   3rd Qu.:31.06   3rd Qu.: 77.0  
##  Max.   :1.0000   Max.   :1.000   Max.   :4.00   Max.   :41.00   Max.   :45.45   Max.   :100.0  
##    windspeed        count         month          daytime    
##  Min.   : 0.0   Min.   :  1   Min.   : 1.00   Min.   :1.00  
##  1st Qu.: 7.0   1st Qu.: 41   1st Qu.: 4.00   1st Qu.:1.00  
##  Median :13.0   Median :144   Median : 7.00   Median :2.00  
##  Mean   :12.8   Mean   :191   Mean   : 6.53   Mean   :2.42  
##  3rd Qu.:17.0   3rd Qu.:283   3rd Qu.:10.00   3rd Qu.:4.00  
##  Max.   :57.0   Max.   :977   Max.   :12.00   Max.   :4.00
summary (bike_test)
##     holiday         workingday      weather         temp           atemp          humidity  
##  Min.   :0.0000   Min.   :0.00   Min.   :1.0   Min.   : 0.82   Min.   : 1.51   Min.   :  0  
##  1st Qu.:0.0000   1st Qu.:0.00   1st Qu.:1.0   1st Qu.:13.94   1st Qu.:16.66   1st Qu.: 47  
##  Median :0.0000   Median :1.00   Median :1.0   Median :20.50   Median :24.24   Median : 62  
##  Mean   :0.0285   Mean   :0.69   Mean   :1.4   Mean   :20.13   Mean   :23.54   Mean   : 62  
##  3rd Qu.:0.0000   3rd Qu.:1.00   3rd Qu.:2.0   3rd Qu.:26.24   3rd Qu.:31.06   3rd Qu.: 78  
##  Max.   :1.0000   Max.   :1.00   Max.   :3.0   Max.   :39.36   Max.   :44.70   Max.   :100  
##    windspeed        count         month          daytime    
##  Min.   : 0.0   Min.   :  1   Min.   : 1.00   Min.   :1.00  
##  1st Qu.: 7.0   1st Qu.: 46   1st Qu.: 3.00   1st Qu.:1.00  
##  Median :13.0   Median :148   Median : 7.00   Median :2.00  
##  Mean   :12.8   Mean   :193   Mean   : 6.48   Mean   :2.43  
##  3rd Qu.:17.0   3rd Qu.:288   3rd Qu.: 9.00   3rd Qu.:4.00  
##  Max.   :57.0   Max.   :897   Max.   :12.00   Max.   :4.00

Вывод: статистический анализ переменных показал, что пустые ячейки отсутствуют, а также существуют возможные выбросы для переменных “windspeed” и “count” в связи большими отклонениями максимальных значений от 3 квартилей. Тем не менее, эти выбросы могут быть объяснены моделью и не являться для неё преградой, поэтому будет проведена отдельная работа с выбросами, которая выявит их значимость для работы моделей.

Также можно заметить, что в датасетах присутствуют 5 качественых переменных (holiday, workingday, weather, month и daytime) и 5 количественных переменных (temp, atemp, humidity, windspeed, count). Качественные переменные выражены количественными категориями, поэтому проведения работ по преобразованию их в “долгий” текстовый тип не нужно - они готовы к работе без дополнительных преобразований.Среди качественных переменных выбросов не было обнаружены, все их количественные предикторы лежат в опредеённых диапозонах, что исключает возможность опечатки при их заполнении.

Диаграммы рессеяния и коэфициенты корреляции

Вывод: На диаграммах рассеяния видна зависимость между всеми параметрами. Коэффициенты корреляции ненулевые и статистически значимы. Наибольшее влияние на целевую переменную оказывают время суток,температура воздуха (по ощущениям тоже) и влажность. Наибольшая корреляция между переменными - temp и atemp, но это очевидно, ведь температура воздуха по ощущениям в наибольшей степени зависит от фактической. Влажность воздуха имеет достаточно большой коэффициент корреляции с типом погоды, что также может быть легко объяснено, в то время как скорость ветра в некоторой степени зависит от влажности, а значит и от погоды. Месяц года имеет хороший коэффициент корреляции с температурой воздуха и его влажностью, что также очевидно, ведь в разные времена года разные темепературные режимы и состояние воздуха в связи с их цикличностью. В остальных случаях, поверхностного анализа недостаточного - необходим более детальный анализ, либо проверка работы моделей с взаимосвязями переменных, что и будет проделано в рамках этого проекта.

3) Ручные регресионные модели:

Модель №1: ‘count’ ~ ‘daytime’

model_1 <- lm(count ~ daytime, data = bike_train)
summary(model_1)
## 
## Call:
## lm(formula = count ~ daytime, data = bike_train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -293.8  -91.0  -51.7   68.3  739.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    27.38       4.10    6.68  2.5e-11 ***
## daytime        67.60       1.53   44.23  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 164 on 8706 degrees of freedom
## Multiple R-squared:  0.183,  Adjusted R-squared:  0.183 
## F-statistic: 1.96e+03 on 1 and 8706 DF,  p-value: <2e-16

Вывод: время суток и число актов аренды показывали наибольшую корреляцию, поэтому было решено построить модель с такими зависимостями первой. Данная модель регрессии статистически значимая, так как p-value = 2e-16 (<0.05).

Проверка параметров модели:

glance(model_1) 
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik     AIC    BIC  deviance df.residual
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>   <dbl>  <dbl>     <dbl>       <int>
## 1     0.183         0.183  164.     1956.       0     2 -56774. 113554. 1.14e5    2.35e8        8706

Вывод: модель смогла объяснить только 18.3% наблюдений, и её стандартная ошибка равна 164, что говорит о том, модель очень плохо справляется своими задачами и нуждается в замене альтернативами.

Модель №2: ‘count’ ~ ‘daytime’ + ‘temp’ + ‘atemp’+ ‘humidity’

t1<- car::powerTransform(count ~ temp + atemp + daytime + humidity, data = bike_train)
coef(t1)
##       Y1 
## 0.273731
model_2 <- lm(count^0.274 ~ temp + atemp +I(atemp^2) + daytime +I(daytime^2) + humidity, data = bike_train)
summary(model_2)
## 
## Call:
## lm(formula = count^0.274 ~ temp + atemp + I(atemp^2) + daytime + 
##     I(daytime^2) + humidity, data = bike_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1104 -0.5598 -0.0489  0.5207  2.5775 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.422615   0.079566  -17.88   <2e-16 ***
## temp          0.017194   0.006738    2.55    0.011 *  
## atemp         0.085848   0.007242   11.85   <2e-16 ***
## I(atemp^2)   -0.001291   0.000114  -11.30   <2e-16 ***
## daytime       3.304186   0.044506   74.24   <2e-16 ***
## I(daytime^2) -0.553026   0.008679  -63.72   <2e-16 ***
## humidity     -0.006951   0.000481  -14.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.774 on 8701 degrees of freedom
## Multiple R-squared:  0.634,  Adjusted R-squared:  0.634 
## F-statistic: 2.52e+03 on 6 and 8701 DF,  p-value: <2e-16

Вывод: Следующая модель - это модель с несколькими переменными, которые имеют наибольшие коэффициенты корреляции с целевой. Все, кроме одного коэффициента в модели, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).

Построение графиков остаточной диагностики:

autoplot (model_2)

Вывод:

  1. Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)

  2. Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.

  3. Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.

  4. Выбросы есть, но в целом модель не объясняет их поведение и сильно подстраивается под них.

Модель 3: “Полная”

Перебирать отдельно набор переменных - долго. Воспользуемся построением полной модели.

В этой модели будет изучено влияние одновремнно всех переменных выборки на целевую.

model_3 <- lm(count ~ ., data = dplyr::select(bike_train, `holiday`:`daytime`))
summary (model_3)
## 
## Call:
## lm(formula = count ~ ., data = dplyr::select(bike_train, holiday:daytime))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -308.1  -96.5  -30.1   55.4  686.3 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.91733    9.83309    0.70     0.48    
## holiday     -3.02912    9.71556   -0.31     0.76    
## workingday   0.00856    3.47119    0.00     1.00    
## weather     -3.81647    2.75132   -1.39     0.17    
## temp         0.48568    1.27233    0.38     0.70    
## atemp        5.98645    1.17103    5.11  3.3e-07 ***
## humidity    -2.19671    0.10113  -21.72  < 2e-16 ***
## windspeed    0.34126    0.21004    1.62     0.10    
## month        7.31759    0.48787   15.00  < 2e-16 ***
## daytime     50.29471    1.44187   34.88  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 146 on 8698 degrees of freedom
## Multiple R-squared:  0.353,  Adjusted R-squared:  0.353 
## F-statistic:  528 on 9 and 8698 DF,  p-value: <2e-16

Вывод: было решено сделать общую модель со всеми влиятельными переменными. Полная модель оказалась значима (p-value: <2e-16), несмотря на то, что некоторые её коэффициенты не значимы, что говорит о целесообразности её дальнейшего тестирования.

glance (model_3)
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik     AIC    BIC  deviance df.residual
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>   <dbl>  <dbl>     <dbl>       <int>
## 1     0.353         0.353  146.      528.       0    10 -55759. 111539. 1.12e5    1.86e8        8698

Вывод: модель опять чуть лучше справилась со своей задачей - она объяснила 35.3% всех наблюдений, и её стандартная ошибка = 146.

Построение графиков остаточной диагностики:

autoplot (model_3)

Вывод:

  1. Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)

  2. Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.

  3. Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.

  4. Выбросы есть, но в целом модель не объясняет их поведение и сильно подстраивается под них.

Перебор с помощью “пустой” модели

Ручной перебор всех 9 переменных с целевой занял бы крайне много времени, поэтому можно оставить эту задачу компьютеру, а потом посмотреть, какие же сочетания давали наилучшие результаты. Воспользуемся информационным критерием AIc, который покажет нам, какие лучшие модели ему удалось выявитЬ.

model_empty <- lm(count ~ 1, data = dplyr::select(bike_train, `holiday`:`daytime`))

m_forward <- stepAIC(model_empty, 
                     scope = list(lower = model_empty, upper = model_3),
                     direction = 'forward')
## Start:  AIC=90602.3
## count ~ 1
## 
##              Df Sum of Sq       RSS   AIC
## + daytime     1  52720548 234641764 88839
## + temp        1  44730893 242631419 89131
## + atemp       1  44296295 243066017 89146
## + humidity    1  28484312 258878000 89695
## + month       1   7590617 279771696 90371
## + weather     1   5171141 282191171 90446
## + windspeed   1   2682806 284679506 90523
## <none>                    287362312 90602
## + workingday  1     18893 287343420 90604
## + holiday     1      6551 287355761 90604
## 
## Step:  AIC=88839.3
## count ~ daytime
## 
##              Df Sum of Sq       RSS   AIC
## + temp        1  32045099 202596665 87563
## + atemp       1  31948909 202692855 87567
## + humidity    1  11973736 222668028 88385
## + month       1   8013594 226628170 88539
## + weather     1   4201890 230439874 88684
## + windspeed   1    313083 234328681 88830
## <none>                    234641764 88839
## + workingday  1     27238 234614526 88840
## + holiday     1      6432 234635332 88841
## 
## Step:  AIC=87562.6
## count ~ daytime + temp
## 
##              Df Sum of Sq       RSS   AIC
## + humidity    1  11127784 191468881 87073
## + weather     1   3076055 199520610 87431
## + month       1   2010581 200586084 87478
## + windspeed   1    678145 201918519 87535
## + atemp       1    174827 202421838 87557
## <none>                    202596665 87563
## + holiday     1      5935 202590730 87564
## + workingday  1       201 202596464 87565
## 
## Step:  AIC=87072.7
## count ~ daytime + temp + humidity
## 
##              Df Sum of Sq       RSS   AIC
## + month       1   5009563 186459319 86844
## + atemp       1    683656 190785226 87044
## + weather     1    161031 191307850 87067
## <none>                    191468881 87073
## + windspeed   1     23090 191445791 87074
## + holiday     1      3553 191465328 87074
## + workingday  1      2411 191466470 87075
## 
## Step:  AIC=86843.8
## count ~ daytime + temp + humidity + month
## 
##              Df Sum of Sq       RSS   AIC
## + atemp       1    527071 185932248 86821
## + weather     1     43899 186415419 86844
## <none>                    186459319 86844
## + holiday     1      5146 186454173 86846
## + windspeed   1      2790 186456528 86846
## + workingday  1       213 186459106 86846
## 
## Step:  AIC=86821.1
## count ~ daytime + temp + humidity + month + atemp
## 
##              Df Sum of Sq       RSS   AIC
## + windspeed   1     45077 185887171 86821
## <none>                    185932248 86821
## + weather     1     29738 185902510 86822
## + holiday     1      1982 185930266 86823
## + workingday  1        23 185932225 86823
## 
## Step:  AIC=86821
## count ~ daytime + temp + humidity + month + atemp + windspeed
## 
##              Df Sum of Sq       RSS   AIC
## <none>                    185887171 86821
## + weather     1     40979 185846192 86821
## + holiday     1      2027 185885144 86823
## + workingday  1        12 185887159 86823

Вывод: наименьшим AIC = 86821 (Информационный критерий Акаике для статистики качества приближения) обладает подобранная модель С переменными: время суток, температура,температура воздуха, влажность, месяц, температрура по ощущениям (+скорость ветра). Стоит проанализировать эту модель для проверки качества её работы.

Модель №4: count ~ daytime + temp + humidity + month + atemp + windspeed

model_4 <- lm(count ~ daytime + temp + humidity + month + atemp + windspeed, data = bike_train)
summary(model_4)
## 
## Call:
## lm(formula = count ~ daytime + temp + humidity + month + atemp + 
##     windspeed, data = bike_train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -316.4  -96.4  -29.9   55.0  689.1 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.6800     9.5002    0.60     0.55    
## daytime      50.1097     1.4355   34.91  < 2e-16 ***
## temp          0.4384     1.2702    0.35     0.73    
## humidity     -2.2591     0.0907  -24.91  < 2e-16 ***
## month         7.3659     0.4865   15.14  < 2e-16 ***
## atemp         6.0360     1.1692    5.16  2.5e-07 ***
## windspeed     0.3024     0.2082    1.45     0.15    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 146 on 8701 degrees of freedom
## Multiple R-squared:  0.353,  Adjusted R-squared:  0.353 
## F-statistic:  792 on 6 and 8701 DF,  p-value: <2e-16

Вывод: все, кроме двух коэффициентов в модели, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).

Проверка параметров модели:

glance (model_4)
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik     AIC    BIC  deviance df.residual
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>   <dbl>  <dbl>     <dbl>       <int>
## 1     0.353         0.353  146.      792.       0     7 -55760. 111535. 1.12e5    1.86e8        8701

Вывод: модель уже лучше справилась со своей задачей - она объяснила по-пренему объяснила 35,3% всех наблюдений, но её стандартная ошибка = 146.

Построение графиков остаточной диагностики:

autoplot (model_4)

Вывод:

  1. Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)

  2. Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.

  3. Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.

  4. Выбросы есть, но в целом модель не объясняет их поведение и сильно подстраивается под них.

Модель №5: count ~ daytime + humidity + month + atemp

В этой модели му удалили незначимую переменную из прошлой модели и проверяем изменение работы модели после этого действия.

model_5 <- lm(count ~ daytime + humidity + month + atemp, data = bike_train)
summary(model_5)
## 
## Call:
## lm(formula = count ~ daytime + humidity + month + atemp, data = bike_train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -319.7  -96.9  -29.9   55.8  691.2 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.8944     8.2387    1.44     0.15    
## daytime      50.2830     1.4306   35.15   <2e-16 ***
## humidity     -2.2988     0.0872  -26.37   <2e-16 ***
## month         7.3071     0.4848   15.07   <2e-16 ***
## atemp         6.4143     0.1941   33.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 146 on 8703 degrees of freedom
## Multiple R-squared:  0.353,  Adjusted R-squared:  0.353 
## F-statistic: 1.19e+03 on 4 and 8703 DF,  p-value: <2e-16

Вывод: все коэффициенты в модели значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).

Проверка параметров модели:

glance (model_5)
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik     AIC    BIC  deviance df.residual
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>   <dbl>  <dbl>     <dbl>       <int>
## 1     0.353         0.353  146.     1187.       0     5 -55761. 111534. 1.12e5    1.86e8        8703

Вывод: модель также справилась со своей задачей - она объяснила по-пренему объяснила 35,3% всех наблюдений, но её стандартная ошибка = 146.

Построение графиков остаточной диагностики:

autoplot (model_5)

Вывод:

  1. Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики)

  2. Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.

  3. Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.

  4. Выбросы есть, и они хуже объясняются, чем в предыдущий модели.

Изучение взаимосвязей между переменными для построения моделей с учётом взаимосвязей факторов

Теперь отойдём от применения простых регрессионных моделей и рассмотрим более сложные. Для визуализации нескольких регресионных прямых, преобразуем рассамтриваемые переменные в отдельных данных в качественные:

bike_train1 <- bike_train 
  bike_train1$daytime <- as.character(bike_train1$daytime) 

bike_train1 <- bike_train1
  bike_train1$weather <- as.character(bike_train1$weather) 

bike_train1 <- bike_train1
      bike_train1$month <- as.character(bike_train1$month) 
        
bike_train1 <- bike_train1
      bike_train1$workingday <- as.character(bike_train1$workingday)

Теперь построим визуализации регресионных прямых:

dots1 <- bike_train1 %>%
  ggplot(aes(count, temp, color = workingday)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title="1", colour = "Рабочие дни", x = "Число актов аренды", y = "Температура воздуха")

dots2 <-bike_train1 %>%

  ggplot(aes(count, temp, color = daytime)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title="2",colour = "Время суток", x = "Число актов аренды", y = "Температура воздуха")

dots3 <-bike_train1 %>%

  ggplot(aes(count, temp, color = month)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title="3", colour = "Месяц", x = "Число актов аренды", y = "Температура воздуха")

dots4 <-bike_train1 %>%

  ggplot(aes(count, temp, color = weather)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = 'lm', se = F)+
  labs(title="4",colour = "Состояние погоды", x = "Число актов аренды", y = "Температура воздуха")

grid.arrange(dots1, dots2, dots3, dots4, ncol=2)

Вывод:

1 график: Какой-либо зависимости между переменными незамечено. Можно сказать, что происходит много числов актов аренды при высокой температуре воздуха в выходные данные, но данная закономерность выражена слаба. Всё же никаких взаимосвязей между этими переменными обнаружено не было;

2 график: Можно заявить, что независимо от температуры воздуха, в ночное время суток берётся в аренду крайне низкое число велосипедов. В остальные времена суток - распределение актов аренды равномерно и не зависит от температуры воздуха;

3 график: Можно заметить, что в зимние времена года, когда присутствовует низкая температура воздуха, число актов аренды минимально. В остальных месяцах число актов аренды распределено равномерно без подчинения к температуре воздуха;

4 график: Сразу бросается в глаза тот факт, что большинство актов аренды происходят во время погоды, при которой ясно, малооблачно или немного облачно. Крайне мало велосипедов арендавалось во время ливня, града, грома, тумана и снега. Другими словами, состояние погоды влияет на число актов аренды по ухудшению ей условий (1 тип -наибольшее число, 4 тип - наименьшее число);

Модель №6: count ~ temp + daytime + humidity + daytimeworkingday + month + weatherdaytime

model_6<- lm(count ~ temp +  daytime + humidity + daytime*workingday + month + weather*daytime, data = bike_train)
summary(model_6)
## 
## Call:
## lm(formula = count ~ temp + daytime + humidity + daytime * workingday + 
##     month + weather * daytime, data = bike_train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -324.6  -97.2  -30.3   55.0  693.4 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          7.5730    12.3121    0.62  0.53851    
## temp                 6.8673     0.2114   32.48  < 2e-16 ***
## daytime             55.8322     3.8841   14.37  < 2e-16 ***
## humidity            -2.1476     0.0967  -22.20  < 2e-16 ***
## workingday         -27.8581     7.8322   -3.56  0.00038 ***
## month                7.3488     0.4860   15.12  < 2e-16 ***
## weather             17.6532     5.7794    3.05  0.00226 ** 
## daytime:workingday  11.2582     2.9147    3.86  0.00011 ***
## daytime:weather     -8.9447     2.1408   -4.18    3e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 146 on 8699 degrees of freedom
## Multiple R-squared:  0.354,  Adjusted R-squared:  0.353 
## F-statistic:  595 on 8 and 8699 DF,  p-value: <2e-16

Вывод: все коэффициенты в модели значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).

Проверка параметров работы моделей:

glance(model_6)
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik     AIC    BIC  deviance df.residual
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>   <dbl>  <dbl>     <dbl>       <int>
## 1     0.354         0.353  146.      595.       0     9 -55756. 111532. 1.12e5    1.86e8        8699

Вывод: модель объяснила 35.4% наблюдений, а стандартное отклонение = 146. Пока что эта модель с взаимосвязями между переменными показала наилучшие результаты.

Диагностические графики:

autoplot (model_6)

Вывод:

  1. Рассеяние точек говорит о наличии нелинейной связи между переменными (линейная связь накладывается на пунктирную прямую в нуле), но тем не менее, коэффициенты нелинейной модели, основываясь на графике, не критически велики).

  2. Это квантильный график для остатков. При справедливости допущения остатки на графике должны быть близкими к наклонной прямой линии. В данном случае, остатки относительно близко прилегают к прямой (особенно в центре), а некоторые наблюдения достаточно несильно отклоняются от нее. Это говорит о корректности построенной регрессионной модели.

  3. Если дисперсия остатков не постоянна, то облако точек будет искривленным или иметь воронкообразную форму. Здесь же график отдалённо напоминает прямую, что значит, что предположение о гомоскедастичности в целом не нарушено.

  4. Левередж = 3*10/8708 = 0.0035 Выбросы есть, но в целом модель несильно подстраивается под них. Достаточно много наблюдений выходят за лопустимые границы, т.е. являются весомыми при построении модели. Возможно, стоит произвести анализ выбросов.

Исследование работ моделей без выбросов и без гетеродоксичности

model_6a <- model_6 %>% 
  fortify() %>% # сохранили остатки, предсказания и показатели влиятельности в набор
  merge(., bike_train) # добавили исходные данные (merge удаляет дублирующиеся при слиянии столбцы)


ggplot(model_6a, aes(.fitted, count)) +
  geom_point(aes(size = .cooksd, colour = .hat)) +
  geom_line(aes(.fitted, .fitted), 
            colour = "red", linetype = "dashed") +
  scale_colour_gradient(low = "black", high = "red") +
  ggrepel::geom_label_repel(
    data = filter(model_6a, .cooksd > 4/nrow(bike_train)),
    mapping = aes(label = "danger"), alpha = 0.25) +
  labs(title = "Влиятельные наблюдения",
       x = "Прогноз count", y = "Факт count", 
       colour = "Левередж", size = "Расстояние Кука")

Модель №7: count ~ temp + daytime + humidity + daytimeworkingday + month + weatherdaytime (без выбросов и гетеродоксичности)

Влиятельные наблюдения помечены словом “danger”, что говорит об опасности использования данных наблюдений в моделе, которые негавтивно на нее влияют. Уберем выбросы и пересчитаем лучшие прогнозы и ошибки моделей:

model_6a$.cooksd <- as.numeric(model_6a$.cooksd)
model_6a$.cooksd[is.na(model_6a$.cooksd)] <- mean(model_6a$.cooksd, na.rm = T)

bike_train_clean <- model_6a %>%
  filter(.cooksd < 4/nrow(model_6a))

Определение степени целевой переменной для устранения гетеродоксичности:

t3<- car::powerTransform(count ~ temp +  daytime + humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
coef(t3)
##       Y1 
## 0.308175

Построение модели:

model_6b<- lm(count^0.308 ~ temp +  daytime + I(daytime^2)+ workingday + I(month^2)+ humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
summary(model_6b)
## 
## Call:
## lm(formula = count^0.308 ~ temp + daytime + I(daytime^2) + workingday + 
##     I(month^2) + humidity + daytime * workingday + month + weather * 
##     daytime, data = bike_train_clean)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.254 -0.669 -0.051  0.603  3.112 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -1.540505   0.102399  -15.04  < 2e-16 ***
## temp                0.035696   0.002358   15.14  < 2e-16 ***
## daytime             4.206833   0.061073   68.88  < 2e-16 ***
## I(daytime^2)       -0.695317   0.010719  -64.87  < 2e-16 ***
## workingday         -0.438254   0.048477   -9.04  < 2e-16 ***
## I(month^2)         -0.011797   0.001643   -7.18  7.7e-13 ***
## humidity           -0.008818   0.000666  -13.23  < 2e-16 ***
## month               0.226562   0.022873    9.91  < 2e-16 ***
## weather             0.029493   0.036046    0.82     0.41    
## daytime:workingday  0.136210   0.018206    7.48  8.1e-14 ***
## daytime:weather    -0.079945   0.013671   -5.85  5.2e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.907 on 8443 degrees of freedom
## Multiple R-squared:  0.684,  Adjusted R-squared:  0.684 
## F-statistic: 1.83e+03 on 10 and 8443 DF,  p-value: <2e-16

Вывод: все коэффициенты в модели, кроме одного, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).

Диагностические графики:

autoplot(model_6b)

Вывод: практически все диагностические графики демонстрируют превосходную работу модели. Самое главное, что нам удалось добиться практически полной горизонтальной прямой на 3-ем графике.

Модель №8: count ~ atemp + daytime + humidity + daytimeworkingday + month + weatherdaytime (без выбросов и без гетеродокчиности)

Поменяем в моделе без выбросов переменную temp на переменную atemp, которая обладала большим коэффициентом корреляции при разведоточном анализе:

t4<- car::powerTransform(count ~  atemp + daytime + humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
coef(t4)
##       Y1 
## 0.308446

Построение модели:

model_6b1<- lm(count^0.308 ~  atemp + daytime+I(daytime^2) + humidity + daytime*workingday + month + weather*daytime, data = bike_train_clean)
summary(model_6b1)
## 
## Call:
## lm(formula = count^0.308 ~ atemp + daytime + I(daytime^2) + humidity + 
##     daytime * workingday + month + weather * daytime, data = bike_train_clean)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.214 -0.667 -0.055  0.606  3.179 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -1.373797   0.098361  -13.97  < 2e-16 ***
## atemp               0.045904   0.001235   37.17  < 2e-16 ***
## daytime             4.125764   0.060305   68.41  < 2e-16 ***
## I(daytime^2)       -0.683799   0.010628  -64.34  < 2e-16 ***
## humidity           -0.008710   0.000662  -13.16  < 2e-16 ***
## workingday         -0.436680   0.048559   -8.99  < 2e-16 ***
## month               0.063237   0.003063   20.65  < 2e-16 ***
## weather             0.009930   0.036011    0.28     0.78    
## daytime:workingday  0.136447   0.018239    7.48  8.1e-14 ***
## daytime:weather    -0.072265   0.013682   -5.28  1.3e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.909 on 8444 degrees of freedom
## Multiple R-squared:  0.683,  Adjusted R-squared:  0.682 
## F-statistic: 2.02e+03 on 9 and 8444 DF,  p-value: <2e-16

Вывод: все коэффициенты в модели, кроме одного, значимые, и сама она статистически значимая p-value < 2e-16 (<0.05).

Модель №9: count ~ atemp + daytime + humidity + daytimeworkingday + monthdaytime + weatherdaytime + windspeed + humiditymonth (без выбросов и без гетеродоксичности)

Определеим степень, в которую необходимо возвести целевую переменную, чтобы сделать нашу лучшую модель гомодоксичной:

t2<- car::powerTransform(count ~  atemp + daytime + humidity + daytime*workingday + month*daytime + weather*daytime + windspeed + humidity*month, data = bike_train_clean)
coef(t2)
##       Y1 
## 0.314036

Построение модели:

model_6b2m<- lm(count^0.314 ~  atemp + daytime + I(daytime^2) + humidity + daytime*workingday + month*daytime + weather*daytime +  humidity*month + humidity*daytime, data = bike_train_clean)
summary(model_6b2m)
## 
## Call:
## lm(formula = count^0.314 ~ atemp + daytime + I(daytime^2) + humidity + 
##     daytime * workingday + month * daytime + weather * daytime + 
##     humidity * month + humidity * daytime, data = bike_train_clean)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.436 -0.698 -0.062  0.630  3.427 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -1.889901   0.146525  -12.90  < 2e-16 ***
## atemp               0.047121   0.001298   36.30  < 2e-16 ***
## daytime             4.349183   0.072610   59.90  < 2e-16 ***
## I(daytime^2)       -0.717278   0.011201  -64.04  < 2e-16 ***
## humidity            0.001762   0.001893    0.93   0.3521    
## workingday         -0.456241   0.050677   -9.00  < 2e-16 ***
## month               0.087614   0.014039    6.24  4.6e-10 ***
## weather            -0.052298   0.040383   -1.30   0.1953    
## daytime:workingday  0.142583   0.019033    7.49  7.5e-14 ***
## daytime:month       0.013635   0.002766    4.93  8.4e-07 ***
## daytime:weather    -0.048311   0.015746   -3.07   0.0022 ** 
## humidity:month     -0.000858   0.000166   -5.16  2.5e-07 ***
## daytime:humidity   -0.002319   0.000591   -3.92  8.8e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.949 on 8441 degrees of freedom
## Multiple R-squared:  0.685,  Adjusted R-squared:  0.684 
## F-statistic: 1.53e+03 on 12 and 8441 DF,  p-value: <2e-16

Вывод: не все коэффициенты в модели значимые, но сама она статистически значимая p-value < 2e-16 (<0.05).

autoplot(model_6b2m)

Вывод: практически все диагностические графики демонстрируют превосходную работу модели. Самое главное, что нам удалось добиться практически полной горизонтальной прямой на 3-ем графике.

4) Регресионные модели, основанные на деревьях решений

Модель №10: обычное дерево решений

Подключение необходимых библиотек для работы с деревьями решений:

library(scales) # Процентный формат для осей ggplot
## 
## Attaching package: 'scales'
## The following object is masked from 'package:memisc':
## 
##     percent
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
library(rpart) # Деревья в R: recursive partitioning
library(party) # Деревья в R: conditional inference
## Loading required package: grid
## Loading required package: mvtnorm
## Loading required package: modeltools
## Loading required package: stats4
## 
## Attaching package: 'modeltools'
## The following objects are masked from 'package:memisc':
## 
##     Lapply, relabel
## The following object is masked from 'package:car':
## 
##     Predict
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## 
## Attaching package: 'strucchange'
## The following object is masked from 'package:stringr':
## 
##     boundary
library(rpart.plot) # Визуализация деревьев
library(partykit) # Визуализация деревьев
## Loading required package: libcoin
## 
## Attaching package: 'partykit'
## The following objects are masked from 'package:party':
## 
##     cforest, ctree, ctree_control, edge_simple, mob, mob_control, node_barplot,
##     node_bivplot, node_boxplot, node_inner, node_surv, node_terminal, varimp
library(C50) # Деревья в R: C5.0
library(mlr)
## Loading required package: ParamHelpers
## 'mlr' is in maintenance mode since July 2019. Future development efforts will go into its
## successor 'mlr3' (<https://mlr3.mlr-org.com>).
## 
## Attaching package: 'mlr'
## The following object is masked from 'package:modelr':
## 
##     resample

Построение автоматически подобронного дерева решений на обучающей выборке:

rtree <- rpart(formula = count~., data = bike_train)
summary(rtree)
## Call:
## rpart(formula = count ~ ., data = bike_train)
##   n= 8708 
## 
##          CP nsplit rel error   xerror      xstd
## 1 0.3108764      0  1.000000 1.000241 0.0195569
## 2 0.0907987      1  0.689124 0.689303 0.0155017
## 3 0.0263655      2  0.598325 0.603862 0.0137338
## 4 0.0213274      3  0.571959 0.574431 0.0131302
## 5 0.0114109      4  0.550632 0.561840 0.0129384
## 6 0.0100000      5  0.539221 0.556046 0.0127668
## 
## Variable importance
##  daytime     temp    atemp    month humidity  weather 
##       53       15       15       10        5        1 
## 
## Node number 1: 8708 observations,    complexity param=0.310876
##   mean=191.136, MSE=32999.8 
##   left son=2 (2523 obs) right son=3 (6185 obs)
##   Primary splits:
##       daytime  < 1.5     to the left,  improve=0.3108760, (0 missing)
##       atemp    < 29.925  to the left,  improve=0.1320830, (0 missing)
##       temp     < 22.55   to the left,  improve=0.1044840, (0 missing)
##       humidity < 66.5    to the right, improve=0.0716794, (0 missing)
##       month    < 3.5     to the left,  improve=0.0548867, (0 missing)
##   Surrogate splits:
##       temp  < 4.51    to the left,  agree=0.712, adj=0.006, (0 split)
##       atemp < 3.41    to the left,  agree=0.711, adj=0.004, (0 split)
## 
## Node number 2: 2523 observations
##   mean=32.5513, MSE=1584.38 
## 
## Node number 3: 6185 observations,    complexity param=0.0907987
##   mean=255.826, MSE=31371.2 
##   left son=6 (2717 obs) right son=7 (3468 obs)
##   Primary splits:
##       temp     < 19.27   to the left,  improve=0.1344740, (0 missing)
##       atemp    < 23.105  to the left,  improve=0.1334120, (0 missing)
##       month    < 3.5     to the left,  improve=0.1077000, (0 missing)
##       humidity < 74.5    to the right, improve=0.0399052, (0 missing)
##       weather  < 2.5     to the right, improve=0.0275818, (0 missing)
##   Surrogate splits:
##       atemp     < 23.105  to the left,  agree=0.998, adj=0.996, (0 split)
##       month     < 3.5     to the left,  agree=0.751, adj=0.434, (0 split)
##       humidity  < 85.5    to the right, agree=0.578, adj=0.040, (0 split)
##       windspeed < 21.0012 to the right, agree=0.567, adj=0.014, (0 split)
## 
## Node number 6: 2717 observations,    complexity param=0.0213274
##   mean=182.446, MSE=19869.5 
##   left son=12 (1665 obs) right son=13 (1052 obs)
##   Primary splits:
##       month   < 9.5     to the left,  improve=0.1135250, (0 missing)
##       atemp   < 14.015  to the left,  improve=0.0824801, (0 missing)
##       temp    < 12.71   to the left,  improve=0.0817189, (0 missing)
##       daytime < 3.5     to the right, improve=0.0431468, (0 missing)
##       weather < 2.5     to the right, improve=0.0259138, (0 missing)
##   Surrogate splits:
##       humidity < 93.5    to the left,  agree=0.617, adj=0.011, (0 split)
##       atemp    < 17.045  to the left,  agree=0.613, adj=0.001, (0 split)
## 
## Node number 7: 3468 observations,    complexity param=0.0263655
##   mean=313.316, MSE=32858.5 
##   left son=14 (875 obs) right son=15 (2593 obs)
##   Primary splits:
##       humidity < 72.5    to the right, improve=0.0664875, (0 missing)
##       atemp    < 29.925  to the left,  improve=0.0459233, (0 missing)
##       weather  < 2.5     to the right, improve=0.0327580, (0 missing)
##       temp     < 24.19   to the left,  improve=0.0226283, (0 missing)
##       daytime  < 2.5     to the left,  improve=0.0173916, (0 missing)
##   Surrogate splits:
##       weather < 2.5     to the right, agree=0.795, adj=0.187, (0 split)
## 
## Node number 12: 1665 observations
##   mean=144.694, MSE=14110.9 
## 
## Node number 13: 1052 observations
##   mean=242.196, MSE=23157.8 
## 
## Node number 14: 875 observations
##   mean=232.854, MSE=24604.7 
## 
## Node number 15: 2593 observations,    complexity param=0.0114109
##   mean=340.467, MSE=32721.8 
##   left son=30 (749 obs) right son=31 (1844 obs)
##   Primary splits:
##       daytime    < 2.5     to the left,  improve=0.0386463, (0 missing)
##       humidity   < 44.5    to the right, improve=0.0263598, (0 missing)
##       temp       < 23.37   to the left,  improve=0.0228858, (0 missing)
##       atemp      < 26.895  to the left,  improve=0.0224512, (0 missing)
##       workingday < 0.5     to the right, improve=0.0157984, (0 missing)
## 
## Node number 30: 749 observations
##   mean=284.67, MSE=21002.9 
## 
## Node number 31: 1844 observations
##   mean=363.131, MSE=35703.6
rpart.plot(rtree)

Вывод: верхнее число над процентом (191, 256 и т.д.) - это среднее число арендованных велосипедов, исходя из условий отбора. Данная модель, по сути, способна прогнозировать только 6 значений count для каждого дня. Модель самостоятельно отобрала, по её мнению, самые важные/значимые переменные - daytime, temp, humidity, month. Примерно такие же переменные были отобраны при ручном моделировании. Значит, данной моделе можно доверять. Попробуем увеличить число ветвей - это позволит повысить число прогнозируемых значений, что в теории должно повысить точность. Тем не менее, позже посмтроим работу этой модели на обучающей выборке.

Модель №11: дерево решений с модифицированной cp

Вычисление среднего значения арендованных велосипедов за весь обучающий период:

mean(bike_train$count)
## [1] 191.136

Создание дополнительного столбца с разделение в числах актов аренды велосипедов по признаку “больше или меньшее, чем среднее”:

bike_train_a <- bike_train %>%
  mutate( count1 = ifelse(count>=191, "much", "less"))

График рассеяния актов аренды велосипедов в разрезе времени суток, температуры и с разделением по признаку “больше или меньшее, чем среднее”:

cmap <- c("much" =  "green", "less" =  "red")

ggplot(bike_train_a, aes(x = daytime, y = temp)) +
  geom_jitter(aes(colour = count1), alpha = 0.5) +
  labs(title = "Сравнение числа актов аренды велосипедов и его среднего значения", colour = "Число аренды в день", x = "Время суток", y = "t воздуха") +
  scale_x_continuous(breaks = 1:10) +
  scale_colour_manual(values = cmap)

Вывод: можно заметить, что ночью практически каждое наблюдение содержит число актов аренды меньшее, чем среднее значение за весь обучащий период, а также низкая температура воздуха, независимо от времени суток негативно влияет на желание людей брать велосипеды и они берут их реже обычного.

Построение модифицированного дерева решений:

rtree_overfit <- rpart(count ~ ., 
                          data = bike_train,
                          control = 
                            rpart.control(
                              minsplit = 2, # минимальное число примеров для разбиения
                              cp = 0.00001, # минимальное относительное улучшение для разбиения
                              maxdepth = 20 # максимальная глубина для листа
                            ))

rpart.plot(rtree_overfit)

Вывод: модель получилась чересчур сложной, она нуждается в упращении.

Посмторим на зависимость относительной ошибки от cp:

plotcp(rtree_overfit)

Вывод: можно увидеть, как сложность модели влияет на ее точность. Попробуем установить параметр сложности (ср) равный 0.0003.

Древо решений с оптимальной cp:

rtree_overfit1 <- prune(rtree_overfit, cp = 0.0003)
rpart.plot(rtree_overfit1)

Вывод: модель стала проще, но тем не меннее, по-прежнему сложна. Также её прогноза на обучающей выборке

Модель №12: полное дерево решений

Посмотрим на полную дерево решений по обучающей выборке:

rctree <- ctree(count ~ ., data = bike_train)
plot(rctree)

Вывод: она оказалась так же сложна. На всякий случай проверим точность её прогноза на обучающей выборке.

5) Автоматичские регрессионные модели:

Модель № 13: линейная регрессия

Подключение необходимых библиотек для работы с автоматическими моделями:

library(tidyverse) # трансформация и визуализация данных
library(mlr) # фреймворк для машинного обучения
library(stringr) # работа со строками
library(forcats) # работа с факторами
library(parallelMap) # распараллеливание задач
library(rpart.plot) # визуализация деревьев rpart

Созданим задачу для автоматической модели: целевая переменная count, данные bike:

bike_task <- makeRegrTask(id = "Аренда велосипедов", 
                         data = bike, target = "count",
                         fixup.data = "warn")
print(bike_task)
## Supervised task: Аренда велосипедов
## Type: regr
## Target: count
## Observations: 10886
## Features:
##    numerics     factors     ordered functionals 
##           9           0           0           0 
## Missings: FALSE
## Has weights: FALSE
## Has blocking: FALSE
## Has coordinates: FALSE

Посмотрим, какие модели существуют для автоматического прогнозирования на основе регресии:

listLearners(obj = "regr") %>%
  dplyr::select(name, class) %>% 
  mutate(name = substr(name, 1, 60)) %>%
  head(10) %>%
  arrange(name)
##                                                            name            class
## 1                            Bayesian Additive Regression Trees regr.bartMachine
## 2                                                 Bayesian CART       regr.bcart
## 3                                     Bayesian Gaussian Process         regr.bgp
## 4  Bayesian Gaussian Process with jumps to the Limiting Linear       regr.bgpllm
## 5                                         Bayesian Linear Model         regr.blm
## 6      Bayesian regularization for feed-forward neural networks        regr.brnn
## 7                               Bayesian Treed Gaussian Process        regr.btgp
## 8  Bayesian Treed Gaussian Process with jumps to the Limiting L     regr.btgpllm
## 9                                   Bayesian Treed Linear Model        regr.btlm
## 10                                            Gradient Boosting         regr.bst

Поиск точного колиества таких моделей

listLearners(obj = "regr") %>% count()
## # A tibble: 1 x 1
##       n
##   <int>
## 1    59

Вывод: их количество - 59.

Прогноз по самым популярным и точным моделям: линейная регрессия:

lrn_lm <- makeLearner("regr.lm", id = "Линейная регрессия")
lrn_lm
## Learner Линейная регрессия from package stats
## Type: regr
## Name: Simple Linear Regression; Short name: lm
## Class: regr.lm
## Properties: numerics,factors,se,weights
## Predict-Type: response
## Hyperparameters:

Взаимосвязь коэффициентов модели и её проверка:

m_lm <- train(lrn_lm, bike_task)

m_lm %>% 
  getLearnerModel() %>%
  summary()
## 
## Call:
## stats::lm(formula = f, data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -314.2  -96.5  -30.7   55.9  685.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.0220     8.7741    1.26    0.209    
## holiday      -4.6855     8.6788   -0.54    0.589    
## workingday   -0.0405     3.1053   -0.01    0.990    
## weather      -3.4616     2.4667   -1.40    0.161    
## temp          2.3462     1.0698    2.19    0.028 *  
## atemp         4.2313     0.9846    4.30  1.7e-05 ***
## humidity     -2.1815     0.0902  -24.19  < 2e-16 ***
## windspeed     0.3420     0.1876    1.82    0.068 .  
## month         7.3987     0.4343   17.03  < 2e-16 ***
## daytime      49.6924     1.2903   38.51  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 146 on 10876 degrees of freedom
## Multiple R-squared:  0.352,  Adjusted R-squared:  0.351 
## F-statistic:  656 on 9 and 10876 DF,  p-value: <2e-16

Вывод: некоторые переменные не значимы, но в целом, модель можно использовать, так как сама она статистичсеки значима?: p-value<0.05.

Разделение данных на обучающий и тестовый период для работ с автоматическими моделями:

# Получаем номера строк для обучающей и тестовой выборки
set.seed(123) # для воспроизводимости
bike_n <- getTaskSize(bike_task)
bike_train_task <- sample(bike_n, size = bike_n * 0.7)
bike_test_task <- base::setdiff(1:bike_n, bike_train_task)

# Строим модель только на обучающей выборке
bike_mod_lm <- train(lrn_lm, task = bike_task, 
                    subset = bike_train_task)

Прогноз с помощью неё на обучающем периоде:

bike_pred_lm <- predict(bike_mod_lm, 
                       task = bike_task, 
                       subset = bike_train_task)

bike_pred_lm
## Prediction: 7620 observations
## predict.type: response
## threshold: 
## time: 0.00
##          id truth response
## 2463   2463   207  298.956
## 2511   2511   420  293.396
## 10419 10419   242  178.309
## 8718   8718   526  402.898
## 2986   2986   313  264.311
## 1842   1842   154  277.009
## ... (#rows: 7620, #cols: 3)

Расчет ошибок модели:

reg_ms <- list(rmse, mape, mae)

performance(bike_pred_lm,
            measures = reg_ms) %>%
  round(2)
##   rmse   mape    mae 
## 144.95   3.16 107.44

Вывод:точность модели плохая, поэтому будут использованы альтернативные автоматические модели для их последующего включение в отбор.

Модель: randomforest - также одна из самых популярных

lrn_lm1 <- makeLearner("regr.randomForest")
lrn_lm1
## Learner regr.randomForest from package randomForest
## Type: regr
## Name: Random Forest; Short name: rf
## Class: regr.randomForest
## Properties: numerics,factors,ordered,se,oobpreds,featimp
## Predict-Type: response
## Hyperparameters:

Присвоение модели выборки:

m_lm1 <- train(lrn_lm1, bike_task)

Просмотр показателей модели (переменные), которые сложно интерпретировать:

m_lm1 %>% 
  getLearnerModel() %>%
  summary()
##                 Length Class  Mode     
## call                4  -none- call     
## type                1  -none- character
## predicted       10886  -none- numeric  
## mse               500  -none- numeric  
## rsq               500  -none- numeric  
## oob.times       10886  -none- numeric  
## importance          9  -none- numeric  
## importanceSD        0  -none- NULL     
## localImportance     0  -none- NULL     
## proximity           0  -none- NULL     
## ntree               1  -none- numeric  
## mtry                1  -none- numeric  
## forest             11  -none- list     
## coefs               0  -none- NULL     
## y               10886  -none- numeric  
## test                0  -none- NULL     
## inbag               0  -none- NULL

Разделение данных на обучающую и тестовую выборки и проноз:

# Получаем номера строк для обучающей и тестовой выборки
set.seed(123) # для воспроизводимости
bike_n1 <- getTaskSize(bike_task)
bike_train_task1 <- sample(bike_n1, size = bike_n1 * 0.7)
bike_test_task1 <- base::setdiff(1:bike_n1, bike_train_task1)

# Строим модель только на обучающей выборке
bike_mod_lm1 <- train(lrn_lm1, task = bike_task, 
                    subset = bike_train_task1)

bike_pred_lm1 <- predict(bike_mod_lm1, 
                       task = bike_task, 
                       subset = bike_train_task1)

bike_pred_lm1
## Prediction: 7620 observations
## predict.type: response
## threshold: 
## time: 0.53
##          id truth response
## 2463   2463   207  202.916
## 2511   2511   420  392.187
## 10419 10419   242  262.529
## 8718   8718   526  472.980
## 2986   2986   313  333.613
## 1842   1842   154  201.492
## ... (#rows: 7620, #cols: 3)

Расчет ошибок модели на обучающем периоде:

reg_ms1 <- list(rmse, mape, mae)

performance(bike_pred_lm1,
            measures = reg_ms) %>%
  round(2)
##  rmse  mape   mae 
## 66.60  0.80 44.57

Вывод: Малые ошибки модели связаны с тем, что randomForest в результате обучения испльзует не одно, а несколько деревьев решений, которые по отдельности показывают не лучшие результаты. Точность получилась неплохой не из-за того, что данный метод для оценки хороший, а за счет чересчур множественных ответвлений и условий при построении дерева. Возможно, что тестовый период окажется не таким точным.

Модель уже лучше, по сравнению с lm, но скорее всего, уступает ручным моделям. Тем не менее, она поучаствует в конкурсе.

*Далее, мы решили посчитать проноз по модели “нечетких множеств”. В течение нашего обучения мы часто сталкивались с данным подходом, который помогал решать множесто логистических задач, связанных с неопределенностью. Здесь не будут представлены chunks с этой моделью, т.к. ее выполнение требует очень долгого времени и в целом нагружает компьютер. Однако, ошибки модели будут представлены: rmse: 132.58 mape: 1.87 mae: 94.95 Результаты ошибок не самые лучшие, значит, модель не будет участвовать в борьбе за лучшую

Модели №14: “bcart”

MLR содержит 59 моделей регресии, считать каждую вручную - нецелесообразно. Попробуем отобрать модели, наиболее подходящие к исследуемым данным. Исходя из характеристик, которыми могут быть описаны модели, выделим 2 значимые: 1) SE - позволяет рассчитывать и прогнозировать стандартные ошибки модели. это позволит выбирать лучшую модель. 2) Numeric - позволяет прогнозировать числовые переменные (count - целевая числовая переменная) В остальном данные не имеют каких-либо особых харакетристик, которые позволили бы отобрать подходящие модели.

listLearners(obj = "regr", properties = c("se", "numerics")) %>%
  dplyr::select(name, class) %>% 
  mutate(name = substr(name, 1, 60)) %>%
  head(50) %>%
  arrange(name)
##                                                            name             class
## 1                                                 Bayesian CART        regr.bcart
## 2                                     Bayesian Gaussian Process          regr.bgp
## 3  Bayesian Gaussian Process with jumps to the Limiting Linear        regr.bgpllm
## 4                                         Bayesian Linear Model          regr.blm
## 5                               Bayesian Treed Gaussian Process         regr.btgp
## 6  Bayesian Treed Gaussian Process with jumps to the Limiting L      regr.btgpllm
## 7                                   Bayesian Treed Linear Model         regr.btlm
## 8                                              Gaussian Process        regr.GPfit
## 9                                            Gaussian Processes      regr.gausspr
## 10                                Generalized Linear Regression          regr.glm
## 11                                                      Kriging           regr.km
## 12                           Local Approximate Gaussian Process         regr.laGP
## 13                                                Random Forest regr.randomForest
## 14                                               Random Forests       regr.ranger
## 15                                           Regression Splines          regr.crs
## 16                                     Simple Linear Regression           regr.lm

Вывод: итого было отобрано 16 моделей, соответствующие заданным параметрам.

Подключение необходимых библиотек для работы с отобранными моделями:

library(tgp)
library(GPfit)
library(DiceKriging)
library(laGP)
library(ranger)
library(crs)
## Registered S3 method overwritten by 'crs':
##   method         from
##   predict.gsl.bs np
## Categorical Regression Splines (version 0.15-31)
## [vignette("crs_faq") provides answers to frequently asked questions]

Вывод: lm и randomforest уже расчитаны. Уберем их.

Создание нескольких моделей:

turnover_learners <- 
  makeLearners(cls = c("bcart", "bgp", "bgpllm", "blm", "btgp", "btgpllm", "btlm", "GPfit","gausspr", "glm", "km", "laGP",  "ranger", "crs" ),
               ids = c("bcart", "bgp", "bgpllm", "blm", "btgp", "btgpllm", "btlm", "GPfit","gausspr", "glm", "km", "laGP",  "ranger", "crs"),
               type = "regr", predict.type = "se")

Создаем задачу - обозначаем целевую переменную:

turnover_task <- makeRegrTask(id = "Аренда велосипедов", 
                         data = bike, target = "count",
                         fixup.data = "warn")

Облегчение процесса расчета прогнозов:

set.seed(123, "L'Ecuyer") # При использовании параллельных вычислений необходимо выбирать алгоритм генерации СЧ "L'Ecuyer"

# Определение количества доступных ядер процессора
num_cores <- parallel::detectCores()

# Запуск параллельных вычислений
parallelStartSocket(num_cores) # для Windows
## Starting parallelization in mode=socket with cpus=8.

Вывод прогнозирование (Данный chunk не выполняется до конца - видимо требует слишком много оперативной памяти, компьютеры не справляются. При уменьшении числа моделей в makelearners ситуация не меняется. Было принято решение - рассчитать модели вручную по отдельности. Ошибки моделей будут представлены ниже.)

turnover_bench <- turnover_learners %>% benchmark(tasks = turnover_task, show.info = FALSE)

“bcart”, “bgp”, “bgpllm”, “blm”, “btgp”, “btgpllm”, “btlm”, “GPfit”,“gausspr”, “glm”, “km”, “laGP”, “ranger”, “crs”

Постановка задачи моделям:

lrn_bcart <- makeLearner("regr.bcart")
lrn_bcart
## Learner regr.bcart from package tgp
## Type: regr
## Name: Bayesian CART; Short name: bcart
## Class: regr.bcart
## Properties: numerics,se,factors
## Predict-Type: response
## Hyperparameters:

Присвоение моделям выборки:

m_bcart1 <- train(lrn_lm1, bike_task)

Показатели модели (переменные), котоыре сложно интерпретировать:

m_bcart1 %>% 
  getLearnerModel() %>%
  summary()
##                 Length Class  Mode     
## call                4  -none- call     
## type                1  -none- character
## predicted       10886  -none- numeric  
## mse               500  -none- numeric  
## rsq               500  -none- numeric  
## oob.times       10886  -none- numeric  
## importance          9  -none- numeric  
## importanceSD        0  -none- NULL     
## localImportance     0  -none- NULL     
## proximity           0  -none- NULL     
## ntree               1  -none- numeric  
## mtry                1  -none- numeric  
## forest             11  -none- list     
## coefs               0  -none- NULL     
## y               10886  -none- numeric  
## test                0  -none- NULL     
## inbag               0  -none- NULL

Разделение данных на обучающую и тестовую выборки и прогноз:

# Получаем номера строк для обучающей и тестовой выборки
set.seed(123) # для воспроизводимости
bike_n1 <- getTaskSize(bike_task)
bike_train_task1 <- sample(bike_n1, size = bike_n1 * 0.7)
bike_test_task1 <- base::setdiff(1:bike_n1, bike_train_task1)

# Строим модель только на обучающей
bike_mod_bcart1 <- train(lrn_bcart, task = bike_task, 
                    subset = bike_train_task1)
## 
## burn in:
## **GROW** @depth 0: [5,0.499943], n=(3619,4001)
## **GROW** @depth 1: [5,0.706953], n=(2530,1340)
## **GROW** @depth 2: [5,0.793103], n=(1576,434)
## **GROW** @depth 2: [6,0.43], n=(299,1965)
## **GROW** @depth 3: [3,0], n=(1172,404)
## **GROW** @depth 2: [2,0], n=(135,299)
## **GROW** @depth 2: [6,0.5], n=(983,1889)
## **PRUNE** @depth 2: [2,0]
## **GROW** @depth 3: [9,0], n=(78,326)
## **GROW** @depth 2: [9,0.333333], n=(1208,448)
## **GROW** @depth 3: [8,0.454545], n=(726,482)
## **GROW** @depth 3: [2,0], n=(174,406)
## **GROW** @depth 5: [4,0.489796], n=(580,1041)
## **GROW** @depth 2: [2,0], n=(428,850)
## **GROW** @depth 5: [6,0.56], n=(180,1126)
## **GROW** @depth 3: [6,0.49], n=(176,252)
## **PRUNE** @depth 3: [6,0.53]
## **GROW** @depth 3: [4,0.285714], n=(14,540)
## **GROW** @depth 3: [8,0.545455], n=(496,206)
## r=1000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; n=(83,289,496,206,441,113,14,540,165,366,1518,812,748,472,1068,289)
## **PRUNE** @depth 5: [3,0]
## **GROW** @depth 4: [9,0], n=(146,106)
## **CPRUNE** @depth 1: var=9, val=0.333333->0, n=(521,875)
## **PRUNE** @depth 3: [6,0.46]
## **GROW** @depth 4: [2,0], n=(449,1069)
## **PRUNE** @depth 3: [2,0]
## **GROW** @depth 3: [5,0.155116], n=(126,668)
## **CPRUNE** @depth 0: var=5, val=0.327529->0.344828, n=(1804,5816)
## **GROW** @depth 5: [5,0.551724], n=(206,226)
## **PRUNE** @depth 5: [5,0.551724]
## **GROW** @depth 4: [9,0.666667], n=(290,76)
## **PRUNE** @depth 4: [9,0.666667]
## **GROW** @depth 4: [5,0.293047], n=(476,410)
## **GROW** @depth 5: [6,0.5], n=(218,258)
## r=2000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; n=(518,323,154,218,258,418,534,446,165,367,366,910,1114,502,1038,289)
## 
## Sampling @ nn=0 pred locs:
## **PRUNE** @depth 5: [6,0.5]
## **GROW** @depth 5: [9,0.333333], n=(79,75)
## **GROW** @depth 5: [3,0], n=(560,554)
## **PRUNE** @depth 5: [3,0]
## **GROW** @depth 4: [5,0.793103], n=(117,48)
## **GROW** @depth 5: [6,0.43], n=(149,327)
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 4: [9,0.333333], n=(62,305)
## **GROW** @depth 4: [2,0], n=(337,832)
## r=1000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=7 n=(576,378,102,52,147,329,572,716,380,123,51,76,330,337,832,854,456,1038,271)
## **PRUNE** @depth 5: [5,0.172414]
## **GROW** @depth 6: [5,0.655172], n=(383,655)
## **PRUNE** @depth 6: [5,0.655172]
## **GROW** @depth 6: [5,0.172414], n=(76,253)
## **GROW** @depth 3: [2,0], n=(266,450)
## **GROW** @depth 5: [9,0.666667], n=(87,29)
## r=2000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,76,253,171,503,266,450,379,87,29,58,76,331,337,832,998,413,937,271)
## **PRUNE** @depth 5: [9,0.666667]
## **GROW** @depth 7: [2,0], n=(23,53)
## **PRUNE** @depth 3: [2,0]
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 3: [9,0.333333], n=(97,282)
## **PRUNE** @depth 7: [5,0.137931]
## **GROW** @depth 5: [6,0.79], n=(696,302)
## **GROW** @depth 6: [8,0.545455], n=(339,357)
## **GROW** @depth 4: [6,0.49], n=(214,289)
## **GROW** @depth 4: [2,0], n=(336,828)
## **PRUNE** @depth 4: [6,0.49]
## r=3000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,103,226,171,503,716,30,377,114,71,83,348,336,828,362,377,328,357,880,257)
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 7: [6,0.83], n=(82,21)
## **GROW** @depth 6: [6,0.69], n=(118,53)
## **GROW** @depth 3: [6,0.48], n=(186,530)
## **PRUNE** @depth 3: [6,0.48]
## **PRUNE** @depth 6: [6,0.7]
## **PRUNE** @depth 7: [6,0.81]
## r=4000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,103,226,171,503,716,30,376,111,74,83,349,1164,606,346,399,272,681,257)
## **GROW** @depth 4: [2,0], n=(337,832)
## **PRUNE** @depth 4: [2,0]
## **GROW** @depth 3: [2,0], n=(266,450)
## **GROW** @depth 6: [5,0.155116], n=(32,139)
## **GROW** @depth 6: [9,0.666667], n=(200,144)
## **PRUNE** @depth 6: [5,0.155116]
## r=5000 d=[0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0] [0]; mh=8 n=(576,378,199,103,226,171,503,266,450,26,327,121,58,77,351,1169,549,481,200,144,293,681,271)
## Grow: 12.95%, Prune: 6.928%, Change: 22.75%, Swap: 8.26%
bike_pred_bcart1 <- predict(bike_mod_bcart1, 
                       task = bike_task, 
                       subset = bike_train_task1)

bike_pred_bcart1
## Prediction: 7620 observations
## predict.type: response
## threshold: 
## time: 84.98
##        id truth response
## 9165 9165    40  38.4029
## 652   652   122 158.7389
## 8110 8110    37  38.4029
## 3494 3494   153 270.1145
## 2278 2278   379 385.9544
## 8484 8484   356 270.1145
## ... (#rows: 7620, #cols: 3)

Расчет ошибок моделей:

reg_ms1 <- list(rmse, mape, mae)

performance(bike_pred_bcart1,
            measures = reg_ms) %>%
  round(2)
##   rmse   mape    mae 
## 128.69   1.53  89.90

Вывод: как оказалось, из-за технических особенностей наших компютеров, наша команда не может рассчитать модели и их ошибки из-за очень длительного процесса вычиления. Спрогнозировав 3 модели видно, что в целом цедевая ошибка Mape гораздо выше, чем у моделей рассчитаных вручную. К сожалению, пришлось ограничится только этими моделями.

6) Выбор лучшей модели

Прогноз на обучающем периоде:

model_2_pred <- bike_train %>%
  add_predictions(model_2, var = "Predicted") %>%
  dplyr::select(count, Predicted)%>%
mutate(Predicted = Predicted^(1/0.274))
 
model_3_pred <- bike_train %>%
  add_predictions(model_3, var = "Predicted") %>%
  dplyr::select(count, Predicted)
  
model_4_pred <- bike_train %>%
  add_predictions(model_4, var = "Predicted") %>%
  dplyr::select(count, Predicted)

model_5_pred <- bike_train %>%
  add_predictions(model_5, var = "Predicted") %>%
  dplyr::select(count, Predicted)

model_6_pred <- bike_train %>%
  add_predictions(model_6, var = "Predicted") %>%
  dplyr::select(count, Predicted) 

model_7_pred <- bike_train %>%
  add_predictions(model_6b, var = "Predicted") %>%
  dplyr::select(count, Predicted) %>%
mutate(Predicted = Predicted^(1/0.308))

model_8_pred <- bike_train %>%
  add_predictions(model_6b1, var = "Predicted") %>%
  dplyr::select(count, Predicted) %>%
mutate(Predicted = Predicted^(1/0.308))

model_9_pred <- bike_train %>%
  add_predictions(model_6b2m, var = "Predicted") %>%
  dplyr::select(count, Predicted) %>%
mutate(Predicted = Predicted^(1/0.314))

model_rtree <- bike_train %>%
  add_predictions(rtree, var = "Predicted") %>%
  dplyr::select(count, Predicted)

model_rtree_cp <- bike_train %>%
  add_predictions(rtree_overfit1, var = "Predicted") %>%
  dplyr::select(count, Predicted)

model_rctree <- bike_train %>%
  add_predictions(rctree, var = "Predicted") %>%
  dplyr::select(count, Predicted)

Расчёт ошибок прогнозов, построенных на основе отобранных ранее ручных и “деревянных” моделей:

errors <- function(actual, predicted) {
  error <- actual - predicted
  abs_error <- abs(actual - predicted)
  percent_error <- abs_error/actual
  
  ME <- mean(error, na.rm = T)
  MAE <- mean(abs_error, na.rm = T)
  MAPE <- mean(percent_error, na.rm = T)
  MSE <- mean(error^2, na.rm = T)
  RMSE <- sqrt(MSE)
  BIAS <- mean(sum(error)/sum(actual))
  
  tibble(ME, MAE, MAPE, MSE, RMSE, BIAS)
}

models_row <- 
  bind_rows("2" = model_2_pred, ##ручная
          "3" = model_3_pred, ##ручная
          "4" = model_4_pred, ##ручная
          "5" = model_5_pred, ##ручная
          "6" = model_6_pred, ##ручная
          "7" = model_7_pred, ##ручная
          "8" = model_8_pred, ##ручная
           "9" = model_9_pred, ##ручная
          "10" = model_rtree, ##дерево решений 
          "11" = model_rtree_cp, ##дерево решений 
          "12" = model_rctree, ##дерево решений
         
          .id = "Number")

models_row %>%
  group_by(Number) %>%
  summarise(errors = list(errors(count, Predicted))) %>%
  unnest() %>%
  arrange(MAPE) %>%
  knitr::kable()
Number ME MAE MAPE MSE RMSE BIAS
9 33.9135 89.0040 1.00187 19247.07 138.734 0.177432
7 34.6628 89.0989 1.00367 19335.14 139.051 0.181352
8 34.6190 89.2499 1.00425 19381.08 139.216 0.181122
2 26.7644 91.4387 1.12310 19418.70 139.351 0.140028
11 0.0000 63.0498 1.30552 7971.81 89.285 0.000000
12 0.0000 82.6160 1.34796 14562.13 120.674 0.000000
10 0.0000 93.9165 1.73647 17794.19 133.395 0.000000
6 0.0000 107.6542 3.00664 21329.52 146.046 0.000000
5 0.0000 107.8052 3.07060 21352.95 146.126 0.000000
4 0.0000 107.7683 3.07074 21346.71 146.105 0.000000
3 0.0000 107.8035 3.07125 21341.75 146.088 0.000000

Рассчёт ошибок автоматических моделей

models_row1 <- 
  bind_rows("randomForest" = performance(bike_pred_lm1,
            measures = reg_ms) ,
            "lm" = performance(bike_pred_lm,
            measures = reg_ms),
            "bcart" = performance(bike_pred_bcart1,
            measures = reg_ms),

        
         
          .id = "Number")

models_row1 %>%
  group_by(Number) %>%
  
  unnest() %>%
  arrange(mape) %>%
  knitr::kable()
Number rmse mape mae
randomForest 66.5993 0.798577 44.5661
bcart 128.6891 1.528816 89.8989
lm 144.9483 3.160999 107.4409

Вывод:

  1. Таким образом, можно выделить “модель-победительницу”на обучающей выборке - это автоматическая модель “randomForest”. Она обладает наименьшей MAPE, MSE, RMSE и MAE по сравнению с остальными отобранными моделями.

  2. Второе место берёт ручная модель №9, которая представляет из себя модифицированную модель №6, проверенной на выборке без значимых выбросов и лишённой гетеродоксичносности. Она показала наименьшую ошибку MAPE c небольшим отрывом от остальных моделей и наименьшие знечения остальных ошибок среди остальных оставшихся моделей.

  3. Третье место забирает модель №7 - это модель №6 с переменой atemp вместо temp. Она обладает чуть большей MAPE, MAE, MSE, RMSE, чем у первых двух моделей в рейтинге, но чуть меньшими ME и BIAS, чем у модели, занявшей второе место.

  4. Остальные модели являются “аутсайдерами”.

7) Прогнозы ручных моделей на тестовой выборке и сравнение ошибок

Построение прогнозов выбранными моделями на тестовой выборке:

model_test_7 <- lm(count^0.308 ~ temp +  daytime + I(daytime^2)+ workingday + I(month^2)+ humidity + daytime*workingday + month + weather*daytime, data = bike_test)
model_test7_pred <- bike_test %>%
  add_predictions(model_test_7, var = "Predicted") %>%
  mutate(Predicted = Predicted^(1/0.308)) %>%
  dplyr::select(count, Predicted)

model_test_9 <- lm(count^0.314 ~  atemp + daytime + I(daytime^2) + humidity + daytime*workingday + month*daytime + weather*daytime +  humidity*month + humidity*daytime, data = bike_test)
model_test9_pred <- bike_test %>%
  add_predictions(model_test_9, var = "Predicted") %>%
  mutate(Predicted = Predicted^(1/0.314)) %>%
  dplyr::select(count, Predicted)

Расчёт ошибок прогнозов моделей №7 и №8 на тестовой выборке:

errors <- function(actual, predicted) {
  error <- actual - predicted
  abs_error <- abs(actual - predicted)
  percent_error <- abs_error/actual
  
  ME <- mean(error, na.rm = T)
  MAE <- mean(abs_error, na.rm = T)
  MAPE <- mean(percent_error, na.rm = T)
  MSE <- mean(error^2, na.rm = T)
  RMSE <- sqrt(MSE)
  BIAS <- mean(sum(error)/sum(actual))
  
  tibble(ME, MAE, MAPE, MSE, RMSE, BIAS)
}

models_row_test <- 
  bind_rows(
          "7" = model_test7_pred, 
          "9" = model_test9_pred,
         
          .id = "Number")

models_row_test %>%
  group_by(Number) %>%
  summarise(errors = list(errors(count, Predicted))) %>%
  unnest() %>%
  arrange(MAPE) %>%
  knitr::kable()
Number ME MAE MAPE MSE RMSE BIAS
7 23.6644 89.3452 1.05954 18234.0 135.033 0.122407
9 22.8688 89.4207 1.07390 18355.1 135.481 0.118291

Прогноз “randomForest” на тестовой выборке:

bike_mod_lm11 <- train(lrn_lm1, task = bike_task, 
                    subset = bike_train_task1)

bike_pred_lm11 <- predict(bike_mod_lm11, 
                       task = bike_task, 
                       subset = bike_test_task1)

reg_ms <- list(rmse, mape, mae)

performance(bike_pred_lm11,
            measures = reg_ms) %>%
  round(2)
##   rmse   mape    mae 
## 111.53   1.31  74.86

Вывод: наилучшей регрессионной моделью на тестовой выброке стала модель №7. Она обладет наименьшей MAPE, совсем незначительно большими значениями других ошибок по сравнению с моделью, занявшей второе место.второе место - модель №9, а третье - автоматическая модель “randomForest”.

Таким образом, удалось определить наилучшую регрессионную модель для прогнозирования, которая показала наилучшие показатели ошибок на тестовой выброке - это модель №9, которая включала в себя множественную регрессию переменных, которая была протестирована на выборке, лишённой значимых выбросов, взаимосвязь переменных, а также избавленной от возможной гетеродокичности. Эта модель заняла третье место по точности прогнозирования на обучающей выборке, но это не помешало ей занять первое место на тестовой выборке. В заключение хотелось бы скзать, что найденная нами регрессионная модель позволяет довольно точно предсказывать значения целевой переменной - число актов аренды велосипедов в зависимости от некоторых внешних факторов, которые определены другими переменными в датасете.